Dielectric response of a polar fluid trapped in a spherical nanocavity 
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We present extensive Molecular Dynamics simulation results for the structure, static and dynam- 
ical response of a droplet of 1000 soft spheres carrying extended dipoles and confined to spherical 
cavities of radii R = 2.5, 3, and 4 nm embedded in a dielectric continuum of permittivity e' > 1. 
The polarisation of the external medium by the charge distribution inside the cavity is accounted 
for by appropriate image charges. We focus on the influence of the external permittivity e' on the 
static and dynamic properties of the confined fluid. The density profile and local orientational order 
parameter of the dipoles turn out to be remarkably insensitive to e' . Permittivity profiles e(r) inside 
the spherical cavity are calculated from a generalised Kirkwood formula. These profiles oscillate in 
phase with the density profiles and go to a "bulk" value et> away from the confining surface; et is 
only weakly dependent on e', except for e' — 1 (vacuum), and is strongly reduced compared to the 
permittivity of a uniform (bulk) fluid under comparable thermodynamic conditions. 

The dynamic relaxation of the total dipole moment of the sample is found to be strongly dependent 
on e' , and to exhibit oscillatory behaviour when e' = 1; the relaxation is an order of magnitude 
faster than in the bulk. The complex frequency-dependent permittivity e(u) is sensitive to e at low 
frequencies, and the zero frequency limit e(uo = 0) is systematically lower than the "bulk" value e;, 
of the static primitivity. 



I. INTRODUCTION 

Following the pioneering work of Debye Q, Kirk- 
wood || , and Onsager , the bulk dielectric response of 
polar materials is by now well understood £| , and dielec- 
tric response is a method of choice for the experimental 
investigation of molecular dynamics in condensed mat- 
ter. Simulations of model polar systems have played a 
key role in our understanding of dielectric fluids [I| , and 
the subtle problems arising in the simulation of finite, 
but periodically repeated samples of such fluids, linked 
to the proper handling of boundary conditions, have been 
clarified in the eig hties | p ,H- However, with few excep- 
tions H, |j| 0, IlllIT^ . Il3| . much less experimental, theo- 
retical and numerical effort has gone into understanding 
the dielectric response of confined polar fluids. 

Consider a fluid of polar molecules trapped in a finite 
or infinite (along one or two directions) cavity surrounded 
by a dielectric material characterised by a permittivity 
different from that of the bulk polar fluid. The global 
dielectric response of the former is determined by the 
fluctuations and relaxation of the overall dipole moment 
of the trapped fluid. The question we wish to address 
is how the dipolar fluctuations are affected by confine- 
ment, i.e. by the presence of a surface separating the 
polar fluid from the dielectric material which surrounds 
the cavity. One can distinguish between two main effects 
due to the presence of such interfaces. The first is purely 
geometric: how are the dipolar fluctuations affected in 
the vicinity of a non-polarisable interface, compared to 
bulk fluctuations? In particular, can one define a mean- 
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ingful local dielectric permittivity tensor *e(r), when the 
polar molecules are restricted to stay on one side of a 
confining surface, with vacuum on the other side. This 
question has been recently addressed in the case of simple 
geometries (slab or spherical cavity) [H, [l2| • The second 
effect arises from the electric boundary conditions which 
must be satisfied when the confining medium is polaris- 
able, and hence characterised by a permittivity e' > 1. In 
this paper we investigate the second effect in the case of 
a simple polar fluid confined to a spherical cavity carved 
out of a dielectric continuum which extends to infinity in 
all directions. Note that in the limit where the system 
consists of a single polar molecule fixed at the centre of 
the spherical cavity, the system reduces to Onsager's cel- 
ebrated model for the calculation of the permittivity of 
a polar material |3|]. 

The model of a polar fluid in a spherical cavity may 
be regarded as a crude representation of dual physical 
situations. One concerns micro-emulsions where inverse 
micelles are nanodroplets of water in oil, which are sta- 
bilised by a monolayer of surfactants. The majority oil 
phase then provides the embedding dielectric medium. 
The conjugate situation is that of a globular macro- 
molecule (i.e. a protein), made up of polar segments, dis- 
solved in water. The connectivity of the macromolecule is 
then crudely accounted for by confining the unconnected 
polar residues to a spherical volume equal to that of the 
cavity. In that case the solvent (water) plays the role of 
the embedding continuum. 



II. MODEL AND SIMULATION 
METHODOLOGY 

Consider a system of N polar molecules confined to 
a spherical cavity of radius R surrounded by an infinite 
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dielectric continuum of permittivity e'. Following related 
earlier work [l^,[T3 | the model which will be investigated 
is one of spherical molecules carrying extended (rather 
than point) dipoles consisting of two opposite charges 
±<7 displaced symmetrically by a distance d/2 from the 
centre of the molecule, such that the absolute dipole mo- 
ment is p = qd. Let fi be the position of the centre of 
the molecule i, and pi be the unit vector along the dipole 
moment of that molecule. The two charges q± = iq are 
than placed at ?i± = ^±5/4,. 

If cj)(r,r') is the electrostatic potential at r 1 due to a 
unit charge at r, taking proper account of the electro- 
static boundary conditions at the surface of the spheri- 
cal cavity, then the total interaction energy of a pair of 
molecules i and j is: 



v(n,rj) = v (\n - fj\) + 



E 

a,0=± 



qaqaWta^^) (l) 



where i>o (r) is the short-range repulsive potential between 
the spherical molecules, which is chosen to be of inverse 
power form as in |T3 |: 



f<j\ n 
v (r) = 4u ( -J 



(2) 



with n = 12 in practice. 

The exact form of 0(r, r") for the spherical geometry 
is derived in the Appendix by solving Poisson's equation 
with the appropriate electrostatic boundary conditions. 
This cumbersome expression is not well adapted to sim- 
ulations, and may be replaced by the approximation: 
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+ (1 - 2k) 



(R/r) 



(3) 



where r* = (R/r) 2 r, e is the permittivity of the empty 
cavity (e = e$ in practice) and k — e'/(e + e'). The 
potential <fi is seen to reduce to the bare Coulomb po- 
tential when e = e', i.e. in the absence of a dielectric 
discontinuity, and reduces to the classic result for a cav- 
ity surrounded by a conductor (metallic boundary con- 
dition e' = 00) 15] . In the approximation J3J the image 
charge is located at the same position as in the metallic 
boundary case, but its weight differs from -1. 

Atoms outside the sphere, making up the dielectric 
continuum of permittivity e', are assumed to interact 
with the dipolar molecules inside the cavity by the short- 
range potential vo(r) in Eq. These atoms are as- 
sumed to be distributed uniformly with a number density 
p, so that the external potential acting on the molecules 
within the cavity is: 



V ex t(r) = P dr' 



J '2 



dd sin 9 



d(j> Au [-) (4) 



gration leads to: 



Vext(r) 



1 

-x 



(n-2)(n-3)(n-4)r 



(n-3)R-r (n-3)R + r 



(5) 



(R-r) n ~ 3 (R + r) n ~ 3 



The external potential goes through a minimum at the 
origin, so that the external force vanishes at r = 0, as 
expected by symmetry. In practice the reduced density 
p* = pa 3 of the dielectric continuum is chosen to be equal 
to 1, thus mimicking a dense medium. 

The coupled classical equations of motion for the trans- 
lations of the molecular centres and the rotations of the 
dipoles were solved by a standard velocity Verlet algo- 
rithm, using the GROMACS Molecular Dynamics (MD) 
package|l6j under constant temperature conditions, im- 
posed by a Berendsen thermostat and with a time-step 
At = 1 fs. The values of the key physical parameters 
are listed in Table [I] Most simulations were carried 
out for samples of N = 1000 molecules and for three 
cavity radii R = 4, 3, and 2.5 nm. Nominal overall 
densities may be estimated as po — 3N/ (47ri?^ ff ) where 
i? e ff < R is an effective radius of the cavity. The latter 
may be estimated from the radial density profiles p(r), 
to be introduced in the following section, by requiring 
p(r = 2R c g — R) = p(r = 0). This makes the effec- 
tive radius of the cavity roughly half a particle diameter 
smaller than the radius at which the dielectric medium 
starts. It is convenient to introduce the following reduced 
variables: 



parameter 






time-step 


At 


1 


fs 


dipole charge 


q 


0.41843035 


c 


charge separation 


d 


0.12190214 


nm 


diameter 


a 


0.36570642 


nm 


energy scale 


u 


1.8476692 


kj/mol 


temperature 


T 


300 


K 


dipole 


P 


2.4500000 


Debye 


mass 


m 


10 


a.m.u. 



TABLE I: Physical parameters as used in the simulations. 
Both charges ±g of the molecule carry a mass of 5 a.m.u. 



Reduced units 




d* 


0.3333 


P* 


2 


r 


0.0278 


T* 


1.35 



where r 2 = r 2 



2rr' cos 6. A straightforward inte- 



TABLE II: Key parameters in reduced units 
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(47reo)M<7 3 



d* 



Pa = Pov 

I 1 



(6) 



i*2 



ma* 



Values of these reduced variables used in the simulations 
are listed in Table [H] Runs extended over several million 
time-steps, corresponding to phase space trajectories of 
several ns. 

The model considered in this paper is not unlike that 
investigated by Senapati and Chandra 8], who used the 
Stockmayer potential for much smaller systems (N ~ 
100), and restricted their calculations to the case k = 0.5, 
i.e. to a cavity surrounded by vacuum. 




FIG. 1: Density profiles p(r) as functions of the radial distance 
r for a system with TV = 1000 and R = 3 nm (p* = 0.53). 
The solid curve is for dipoles (p* = 2) with k — 0.5. The 
different values of k coincide on this scale. The dotted curve 
is the reference for particles without dipoles (p* = 0). 



III. STATIC PROPERTIES 

This section focuses on the results of our MD simu- 
lations for the structure and static dielectric response 
of the model defined in Sect. [D] for a polar fluid in a 
spherical cavity. We have considered embedding dielec- 
tric continua of permittivities e' = 1 (vacuum), 4, 9, 
and co (metal) corresponding to values of the parame- 
ter k = e'/(e + e') = 0.5, 0.8, 0.9, and 1. The structure of 
the polar fluid is conveniently characterised by the radial 
density profile p(r), where r is the distance of the centre 
of a polar molecule from the centre of the cavity; clearly 
p(r) = for r > R. The computed profiles integrate up 
to the total number of polar molecules in the cavity: 



-in 



drp(r)r 2 = N 



(7) 



Profiles obtained for a cavity of radius R = 3 nm, 
N = 1000, p* = 2 and four values of k are com- 
pared in Fig. ^ to the profile corresponding to non-polar 
molecules (p* = 0) under otherwise identical conditions. 
All profiles exhibit the expected layering near the confin- 
ing spherical surface 8]. As expected the layering effect 
is even more pronounced for the smaller cavity radius 
R = 2.5 nm which we also explored (data not shown). 
There are two striking results: the profiles observed for 
the four different values of k are nearly indistinguish- 
able, i.e. the radial structure turns out to be practi- 
cally independent of the polarisability of the confining 
continuum. However the profile p(r) changes substan- 
tially when p* is set equal to zero, i.e. in the absence 
of any electrostatic coupling between molecules and with 
the embedding medium: the layering is seen in Fig. ^ to 
be considerably enhanced, and to extend deeper towards 
the centre of the cavity. These findings are qualita- 
tively confirmed in the cases of the larger (i? = 4 nm) 




FIG. 2: The order parameter ■ T)) as a function of the 

distance to the centre of the cavity (N = 1000, R = 2.5 
nm). The dotted curve shows the corresponding density pro- 
file p*(r) for k — 0.5. 



and smaller (R = 2.5 nm) cavities. The conclusion to 
be drawn here is that the dipolar interactions between 
molecules tend to smooth out the layering imposed by 
the short-range, excluded volume effects. The orienta- 
tions of individual dipoles relative to the normal to the 
surface are characterised by the local order parameters 
(Pi (p ■ f)) r and (P2(p ■ r)) r , where Pi denotes the Ith 
order Legendre polynomial, r and u are the unit vectors 
along the radial vector r and the dipole moment p, and 
the statistical average is taken over dipole configurations 
within spherical shells of radius r and width 5. Because 
of the jli — > — pi inversion symmetry, {P\(p • r)) r is iden- 
tically zero, while {P^ip • f)) r 1S plotted in Fig. [3 as a 
function of r for several values of n. (i"2(A ' ^))r is seen 
to depend very little on k, except very close to the outer 
surface. The order parameter oscillates somewhat out of 
phase with the oscillations in the density profile shown 
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FIG. 3: The average absolute value {\M(r)\) of the total 
dipole moment of all dipoles within a distance r from the 
centre of the cavity, for the various values of k. From top to 
bottom cavities of size R — 4, 3, and 2.5 nm. 



in a frame of the same figure. Near the maxima of the 
latter, which correspond to well defined shells of polar 
molecules, the order parameter is predominantly nega- 
tive, signalling a preferential orientation of the dipoles 
orthogonal to the radial vector r, i.e. the dipoles orient 
preferentially parallel to the confining surface, irrespec- 
tive of the embedding medium, suggesting a vortex-like 
pattern of the confined dipoles. 

Qualitatively similar behaviour is observed for the 
larger cavities (lower overall densities), except that, as 
expected, the oscillations in (P2(A ' ^))r a re less pro- 
nounced. The preferential alignment of dipoles paral- 
lel to the confining surface was also observed in earlier 
work jl7j . 

Consider next the total dipole moment M (r) within a 
sphere of radius r < R. While (M(r)) vanishes again by 
symmetry, the statistical average of the absolute value of 
M (r) shows an interesting behaviour, illustrated in Fig. El 
for the three different pore radii R under investigation. 
For the two lower densities (\M(r)\) is seen to increase 
roughly as N^ 2 i.e. r 3 / 2 , up to r ~ R/2, as one would 
expect if the N r dipole moments within a sphere of radius 
r were uncorrelated. This part of the curve is essentially 
independent of k. Beyond r ~ R/2 the four curves di- 
verge and show some structure for the smaller cavities 
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FIG. 4: The correlation (/t • M) of the individual dipoles with 
the direction of the total dipole moment of all particles in the 
cavity as a function of the distance r from the centre for a 
R = 4 nm cavity, and 0.5 < « < 1. 



(R = 3 and 2.5 nm). For a given r, the value of (|M(r)|) 
increases with increasing k, reflecting an enhanced effect 
of the image dipoles as e' increases (See. Eq. EJ. Al- 
though not evident in the density profiles of Fig. ^ the 
effect of the polarisation of the confining medium is seen 
to have a very significant effect on the dipolar proper- 
ties of the confined system. This is also evident in Fig. 01 
where we plot the average of the projection of the dipole 
moments of the individual dipoles within a shell of radius 
r and small thickness along the total dipole moment M 
of the sample, (/t • M). The projection is small, but pos- 
itive signalling that the individual dipoles tend to align 
along the overall dipole moment. Moreover, it is quasi- 
independent of the distance r, and increases with k, i.e. 
as one moves from vacuum outside the cavity, to a metal- 
lic confining medium. 

We finally turn to the static dielectric permittivity pro- 
file of the confined fluid. This can be related to the 
dipolar fluctuations by linear response theory 12], or by 
measuring the total polarisation induced in the sample 
by the "external" electric field due to a charge placed at 
the origin of the cavity. Let m(r) denote the microscopic 
polarisation density: 



JY 



The overall dipole moment of the sample is then: 



M 



dfrh(r) 



(8) 



(9) 



•Doavit; 



where the integration is over the whole volume of the 
cavity. The linear response result for the permittivity 
profile e(r) is given by the following generalisation °f 
Kirkwood's classical results for the bulk 0: 



(e(r)-l)(2 £ ' + l) 47T/? 
2e' + e(r) ~ 3 



(m(r) • M) 



(m(r)) ■ (M) 
(10) 
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Far from the confining surface bulk behaviour may be 
expected, and replacing fh by M/V, one recovers Kirk- 
wood's formula. Near the surface rotational invariance is 
broken and the permittivity becomes a tensor with lon- 
gitudinal (i.e. parallel to r) and transverse components. 

Note that equation l|10(l is exact provided that there 
exists a local relationship between the polarisation and 
the internal (Maxwell) electric field. It defines a per- 
mittivity profile e(r) which may be expected to go to a 
constant "bulk" value far from the confining surface, as 
will be confirmed by our simulation data, at least at low 
or moderate densities. An approximate method for esti- 
mating such "bulk" values within cavities has been put 
forward by Berendsen [Hi, IHSI ] , but the limitations of this 
method have been illustrated in Ref . [T^| . 

In MD simulations, the correlation function on the 
right hand side of the Eq. fTTijl is estimated by averaging 
over all dipole moments of particles within a spherical 
shell of radius r and width ~ a. In the "external field" 
method, an additional particle is placed at the origin, 
with the (extended) dipole replaced by a simple proton 
charge e at its centre. If one assumes a local relationship 
between the radial polarisation density P(r) — (f ■ m(r)) 
and the local radial electric field E(r) — f ■ E(r) of the 
form: 



P(r) = e oX {r)E(r) = [e(r) - l]E(r) 



(11) 



then, e(r) follows from elementary electrostatics [TJ. Let 
Q(r) = e + Qind(r) be the total charge contained inside 
a sphere of radius r, which is easily estimated from the 
MD simulations for the extended dipole model. As a 
consequence of the divergence theorem, Q-mdir) inside the 
sphere of radius r is related to the polarisation density 
by: 



P(r) = - 



Qind(r) 



while the electric field E{r) is related to Q{r) by: 



(12) 



(13) 




FIG. 5: The radial density profiles p(r) in the presence of a 
charge at the centre of the cavity for k = 0.5 (R = 3 nm, 
fi* —2). The different values of kappa coincide on this scale. 
The dotted curve is the reference profile for particles with- 
out dipole. In the lower part of the figure the corresponding 
charge distribution Qmd(r)/e is shown. 




FIG. 6: e(r) versus the distance r obtained from the fluctua- 
tion formula il l 01 (open symbols) and from the response to a 
central charge 1141 1 (filled symbols) for fi* = 2, R = 4nm, and 
0.5 < k < 1. 



Substitution of and (JT3J in Eq. JTTJ leads to the 
desired estimate: 



e(r) 



1 



e - 47rr 2 P(r) 1 



i(r)/e 



(14) 



The presence of the central charge introduces some dis- 
tortion of the density profiles near the centre of the cav- 
ity, as illustrated in Fig. [SJ An excluded volume zone 
and subsequent layering now appear at small r, but the 
profiles are virtually unchanged for r > R/2 relative to 
the case without central charge. Note that adding the 
additional particle changes the overall density inside the 
cavity by only one part in 1000, so that a comparison 
between fluctuation and response results remains mean- 
ingful. 



Plots of the induced charge Qi n d(r) inside a sphere of 
radius r are shown as a function of r in Fig. for the 
cavity with radius R = 3 nm. Qind(r) "overscreens" the 
external charge e at the centre, at short distances, before 
oscillating around a negative value and going to zero as 
r — > R . At the lower overall density (R = 4 nm), Q m d(r) 
stabilises around a "bulk" value above — e at intermediate 
distances (data not shown) , while no such "bulk" regime 
is observed at higher density (R = 3 nm) . Remarkably 
Q(r) is practically independent of e' (or k). 

The values of e(r) derived from Eq. (|14|) are com- 
pared to the corresponding values estimated from the 
fluctuation formula (|10fl in Fig. for the lower density 
(R = 4 nm), /i* = 2 and four values of k. Clearly 
Eq. I|14fl can only yield physically acceptable results as 
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r [nm] 



FIG. 7: e(r) versus the distance r obtained from the fluctua- 
tion formula (I10P for fi* — 2, R = 3 nm, and 0.5 < K < 1. 




r [nm] 



FIG. 8: e(r) versus the distance r obtained from the fluctu- 
ation formula llOt for fi* = 2, R = 2.5 nm, and re = 0.8, 0.9, 
and 1. 



long as QindM/e > — 1. At small and large r the oscilla- 
tions in Q(r) lead to unphysical values of e(r), signalling 
the break-down of the local assumption Ijlll) , as already 
noted in Ref. 0] f° r the special case k — 0.5 (e' = 1). 
In that case the data shown in Fig. H3 point, however, to 
two surprising findings. First of all e(r) profiles derived 
from Eq. (|14fl appear to be independent of k, due to the 
quasi-independence of Qi n d(r) on k in Fig. [3J On the 
other hand the e(r) values derived from the fluctuation 
formula 1)11) [) for e' > 1 (k = 0.8, 0.9, and 1) are prac- 
tically independent of k, but lie significantly below the 
results for k = 0.5. In other words, the polarisation of 
the surrounding medium leads to a reduction of the per- 
mittivity of the sample inside the cavity. However the 
response to a central charge e appears to be significantly 
non-linear, but independent of k. It would be very diffi- 
cult to measure the response to a smaller external charge 
at the cavity centre, due to an unfavourable signal to 
noise ratio. 

Results for the smaller cavity (R = 3 nm), i.e. higher 



R 


P* 


K 




e(<j = 0) 


4 nm 


0.23 


0.5 


7.6 


4.8 






0.8 


5.7 


4.5 






0.9 


5.6 


4.6 






1.0 


5.7 


4.8 


3 nm 


0.53 


0.5 


19.7 


10.5 






0.8 


11.3 


9.0 






0.9 


11.3 


9.5 






1.0 


12.3 


10.7 


2.5 nm 


0.92 


0.5 




18.2 






0.8 


20.8 


13.6 






0.9 


20.8 


15.1 






1.0 


21.8 


18.4 



TABLE III: The estimated "bulk" permittivity tb and zero- 
frequency permittivity e(o) = 0) values for the different cavity 
sizes R, reduced densities p* , and re. 



overall density are shown in Fig. As is already clear 
from Fig. Eq. p4(l is no longer applicable, because 
the strong oscillations in Qi n d(r) go repeatedly through 
the value — e and a "bulk" -like regime is never reached. 
Hence only the results from the fluctuation formula 1)10(1 
are shown. The behaviour as a function of k qualita- 
tively confirms one of the observations already made at 
the lower density (R = 4 nm) shown in FigEI namely 
that the "bulk" values of e agree within statistical errors 
for e' > 1 (k = 0.8, 0.9, and 1), but are roughly a factor 
2 lower than the value measured for e' = 1. The error 
bars on the latter are much larger than those associated 
with the data for the polarisable embedding medium. 

At the highest density (R = 2.5 nm), e(r) oscillates 
roughly in phase with the density oscillations. A proper 
"bulk" regime is never reached for the N = 1000 particle 
system, but one can extract a rough value of eb around 
which e(r) oscillates. These values, given in Table IIIII 
depend only weakly on k, except for k = 0.5 (vacuum 
outside the cavity), when eb is roughly a factor of two 
larger than for n > 0.5; the very noisy permittivity pro- 
file for k — 0.5 is not shown in Fig. [S] The best esti- 
mates of eb as a function of cavity radius R and of n are 
listed in Table ITTll The permittivity of the confined fluid 
is strongly reduced compared to its value in a uniform 
(bulk) fluid at the same density and temperature. 



IV. RELAXATION 

The dielectric response of a polar sample is charac- 
terised by the frequency-dependent complex dielectric 
permittivity e(u>) = ei(w) + iei(yj). Within the linear 
response regime the latter is determined by the Laplace 
transform of the dynamical response function $mmW 
which relates the induced total dipole moment of the 
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0.4 0.6 

tips] 

FIG. 9: The total dipole autocorrelation function Cmm (£) 
versus time for 7? = 4 nm, /i* = 2, 0.5 < k < 1. 




0.4 0.6 
t[ps] 

FIG. 11: The total dipole autocorrelation function CMMit) 
versus time for R = 2.5 nm, /x* = 2, 0.5 < « < 1. 




0.4 0.6 
t[ps] 

FIG. 10: The total dipole autocorrelation function Cmm (£) 
versus time for 7? = 3 nm, /i* = 2, 0.5 < /t < 1. 



sample to a time-dependent external field 0, H3 : 

(AM(t))= / $MM(t-t')^ext(*')^' (15) 
J —oo 

where $a/m is a scalar for a spherical sample. According 
to the standard rules of linear response |2(j : 

$MM(t) = -P(ti(t) ■ M(0)) 



-f3C M M(t)(M 2 ) 



(16) 



where the dot denotes a time derivative and CMhi{t) is 
the normalised total dipole moment correlation function 
of the unperturbed sample: 

C M„( t )= (jg( ty o)> (in 



The complex susceptibility is: 



Xmm{z)= / <$>MM(ty zt dt (18) 
Jo 



where z = w + is, and: 



lim xmm(z) = XiM + %M 



(19) 



According to the fluctuation-dissipation theorem, a di- 
rect consequence of Eq. (|16|) : 

Xmm(z) = P{M 2 ) [c M u(t = 0) + izCmm{z)] (20) 

This implies that the spectrum of the correlation function 
Cmm if) is related to the imaginary part of the suscepti- 
bility 



Cmm{oj) — — 
Zn 



'CMM^dt 



7TL0 (M 2 ) 



(21) 



Finally, for the system under consideration, i.e. a po- 
lar fluid confined to a spherical cavity surrounded by a 
dielectric continuum of permittivity e', the frequency- 
dependent permittivity of the sample is related to the 
complex susceptibility by [2l| : 



e(w) - 1 2e' + e 



k B T 



e-1 2e' + e(oj) (M 2 ) 



{M 2 



Xmm(u) 
[XiH + 3X 2 M] 



(22) 



where e = e(w = 0) is the static permittivity of the sam- 
ple. 

The static permittivity e is given by the lo — » limit 
of Eq. 122|) which results in: 



(e-l)(2e' + l) 4tt/3(M 2 } 



2e' + e 



(23) 



Thus e(w = 0) is determined by the fluctuation of the 
total dipole moment M of the spherical sample. The 
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FIG. 12: The dipole autocorrelation lunction Cmm(t, t) ver- 
sus time for R = 3 nm, fi* = 2, k = 1 and several values of 




[ps| 



FIG. 13: The dipole autocorrelation function C'MM{r,t) ver- 
sus time for R = 3 nm, p* — 2, k — 0.5 and several values of 
r. 



result differs from the "bulk" value eb determined as ex- 
plained in Sect. Hill from the profile e(r) calculated from 
Eq. (|10[1 . It characterises the global response of the polar 
fluid trapped in the cavity rather than the local response, 
away from the confining surface. The values of e are seen 
to lie systematically below those of £{,, particularly so in 
the case k — 0.5. 

Hence all information required to compute e(u>) is con- 
tained in the total dipole autocorrelation function (|17|l . 
MD results for Cmm (t) are shown in Figures - ^2 for 
the three cavity radii R = 4, 3, and 2.5 nm, and four 
values of k. The correlation functions are seen to relax 
to zero over a time scale of about 1 ps, which is an order 
of magnitude shorter than the relaxation time observed 
in the bulk for the same model ^4|. The most striking 
feature is the strong sensitivity of Cmm{^} to the po- 
larisability of the confining medium, i.e. to k. For all 
three radii, the relaxation is slowest for k = 1 (metallic 
boundary), and becomes faster as k decreases. Marked 
oscillations appear when k = 0.5 (e' = 1) particularly so 
at the highest density (i? = 2.5 nm). Simulations carried 
out on smaller samples of N — 250 dipoles show that 
the relaxation patterns appear to be independent of the 
sample size characterised by R and N, provided the re- 
duced overall density p* = pa 3 is the same. The decrease 
of the relaxation time of Cmm (t) with k agrees with the 
behaviour predicted for a Debye dielectric 0. There is 
no obvious explanation for the oscillation observed for 
k = 0.5. These oscillations are indicative of collective 
behaviour reminiscent of the dipolaron mode observed 
in MD simulations of longitudinal dipolar fluctuations at 
finite wavenumber in bulk model polar fluids |22|.In an 
effort to gain a better understanding of the dipolar relax- 
ation, we have also computed the normalised correlation 
functions: 



Cmm {r, t) 



(M(r,t) -M(r,0)) 
WW) 



(24) 



where M(r,t) is the instantaneous total dipole moment 
of all the molecules contained inside a sphere of radius 
r < R. The MD data are shown in Figures IT21 and IT31 for 
k = 1 and k = 0.5 respectively. In the former (metal- 
lic boundary) case Cmm {f, t) changes only moderately 
with r. This is not totally unexpected, since the "bulk" 
permittivity of the sample is relatively large under these 
conditions (e ~ 12), as seen from Fig. [7| and hence the 
dielectric discontinuity is not too strong relative to the 
metallic embedding medium. The situation is very dif- 
ferent when k — 0.5 (Fig. H3|) . In this case CMM(f,t) 
changes relatively little with r, and decays monotonically 
except when r = R (corresponding to the autocorrelation 
function of the total dipole moment of the sample), when 
Cmm decays much faster and oscillates. Thus the dy- 
namics of the molecular dipole moments inside the outer 
shell, in direct contact with the surface separating the 
confined fluid from vacuum, has a dramatic effect on the 
total dipole correlation function. 

Real and imaginary parts of the complex susceptibility 
(|T§)) are plotted in Fig. ^] for the cavity of radius R = 3 
nm, and 4 values of k. Xi(w) and X2(w) vary with lu in 
a manner reminiscent of bulk behaviour |22j. However 
they are rather sensitive to k, i.e. to the permittivity of 
the surrounding medium. On the contrary the real and 
imaginary parts of the frequency-dependent permittivity 
denned by Eq. (12211 . are remarkably insensitive to k, ex- 
cept for ex(u>) in the static (to — > 0) limit, as illustrated 
in Fig. ^3 This insensitivity to k reflects the fact that, 
contrary to xmm(w), e(uj) measures the response of the 
polar fluid to the local, internal field. 

The resulting Cole-Cole plots, shown in Fig. El differ 
strongly from the semi-circular shape of the simple Debye 
theory, as one might expect. The high-frequency part 
is very insensitive to K, while the differences in the low- 
frequency range reflect the significant differences in static 
values of e. 
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FIG. 14: The real and imaginary part of the complex suscep- FIG. 15: The real and imaginary part of the complex permit- 
tibility as a function of the frequency u (R = 3nm). tivity as a f unction of the frequency w (R= 3nm). 



V. CONCLUSION 

We have reported the first systematic attempt to in- 
vestigate the dependence of the structure, static and dy- 
namic correlations and response of a drop of polar fluid 
confined to a spherical cavity, on the permittivity e' of 
the embedding medium. The work extends earlier in- 
vestigations which were restricted to the case of a non- 
polarisable external medium (e' = 1) 0, by treat- 
ing the interactions between the charge distribution as- 
sociated with the extended dipoles of the trapped fluid, 
and the image charges in an approximate, but accurate 
way, which provides an efficient alternative to the more 
cumbersome variational method proposed elsewhere |23| . 
The present treatment of electrostatic boundary condi- 
tions for the electric field of individual charges within 
the confined sample does not rely on macroscopic reac- 
tion field considerations, but is restricted to the spherical 
geometry. 

The MD simulations were run for samples of N = 1000 
polar molecules confined to cavities of radii R = 4, 3, and 
2.5 nm, which amount to effective densities of p* = 0.23, 
0.53, and 0.92. Some test runs were carried out for a 
smaller sample of N = 250 molecules, and no significant 
iV-dependence was observed. The larger TV = 1000 par- 
ticle system allows a "bulk" regime to be reached within 
a substantial fraction of the accessible volume (say up 
to r ~ R/2), except for the smallest cavity (i.e. high- 
est effective density p* = 0.92). All calculations were 
made with a reduced extended dipole moment pf = 2, 
comparable to that of water. 

The main conclusions to be drawn from our MD data 
may be summarised as follows: 

a) The structural properties embodied in the den- 
sity profiles p(r) and the order parameter profiles 
{P2{p • r)) are remarkably insensitive to the em- 
bedding medium, i.e. to k. The density profiles 
show significantly less structure than their p = 



i 1 1 1 1 1 1 1 1 1 r 



6- 




FIG. 16: Cole-Cole plot: imaginary versus real part of the 
permittivity (R = 3nm) . 



counterparts. The order parameter profiles point 
to a preferential alignment of the dipoles parallel 
to the confining surface, as already reported in ear- 
lier studies of related systems . 

b) The static permittivity profiles e(r) may be cal- 
culated from the generalised Kirkwood fluctuation 
relation 03 12 . or by measuring the polarisation 
profile P(r), or charge profile Q(r) induced by an 
"external" charge placed at the centre of the spher- 
ical cavity (cf. Eq. (|14f> ). The latter method can 
only be implemented at the lowest density (R = 4 
nm) and points to a significant non-linearity com- 
pared to the predictions of the fluctuation for- 
mula, when the "external" charge is the proton 
charge. The fluctuation formula yields oscillatory 
e(r) profiles which essentially reflect the oscillations 
in the density profiles p(r). At the two lower den- 
sities (R — 4 and 3 nm) the oscillations are suf- 
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ficiently damped away from the confining surface 
for a "bulk" regime to be reached inside the cavity 
allowing the definition of a "bulk" permittivity eb- 
The latter turns out to be relatively insensitive to k, 
except for k — 0.5 (cavity surrounded by vacuum), 
which leads to substantially larger values of ef,. At 
the highest density (R — 2.5 nm), the oscillations 
of e(r) extend up to the centre , so that no proper 
"bulk" regime is reached . A rough estimate of ef, 
may be extracted by averaging over oscillations. In 
all cases the "bulk" permittivity inside the cavity 
is strongly reduced relative to the values expected 
for a uniform (genuinely bulk) fluid under compara- 
ble conditions, an observation also made in earlier 
workHIlll. 

c) While the static properties of a confined drop of 
polar fluid turn out the be surprisingly insensi- 
tive to the permittivity of the external medium ex- 
cept when e' — ► 1, the dynamical properties de- 
pend much more on e' (or equivalently k). The 
most striking illustration is the correlation function 
CmmW of the total dipole moment of the sample, 
which relaxes faster when k ~ 0.5 (cf. Figs.l§l- lll[) . 
As also noted in earlier work [l^ on confined water, 
the relaxation of Cmm {t) is much faster than under 
comparable bulk conditions 0] . While such a be- 
haviour may be rationalised by the absence of long- 
range dipolar interactions with distant molecules in 
the case of the confined system, which may lead to 
a lower collective "inertia" compared to the bulk, 
a detailed theoretical interpretation of the observa- 
tion is still lacking. 

d) The complex, dynamical permittivity was esti- 
mated from Eq. (|22|l The corresponding Cole- 
Cole plots differ considerably from the semi-circular 
shape associated with exponential Debye relax- 
ation, as was to be expected from the complex re- 
laxation pattern of the Cmm(£) correlation func- 
tion. While the high frequency regions of the real 
and imaginary parts of the permittivity are quite 
insensitive to the value of k, significant deviations 
occur in the low- frequency regime (cf. Figs. ITU and 
11511 which are reflected in the right-hand parts of 
the Cole-Cole plots in Fig. ^| The zero-frequency 
limits e = e(cu = 0) of the dynamical permittiv- 
ity differ substantially from the "bulk" limits ej of 
the static permittivity profiles e(r), as shown in Ta- 
ble [^TJ which is not surprising since e and ef, mea- 
sure "global" and "local" responses respectively. 

A theoretical analysis of the dynamical properties of con- 
fined polar fluids is left for future work. We also plan to 
explore the static and dynamic properties of polar flu- 
ids in narrow pores (one-dimensional confinement) and 
in slits (two-dimensional confinement). 
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APPENDIX A 

In this appendix we derive the electrostatic potential 
inside a spherical cavity of radius R and permittivity e, 
which results from the polarisation of an infinite medium 
with permittivity e' surrounding the cavity due to an 
internal charge distribution. 

The electrostatic potential at r, due to a single charge 
Q at position d inside the cavity, can be formally ex- 
panded in Legendre Polynomials Pi as: 



Q 1 
B, 



J2Air l Pi{cos9); r<R 



1=0 



(Al) 



$out -E^( cos 



r>R 



2=0 



where cos — f ■ d, and the expansion coefficients Ai and 
E>\ follow from the boundary conditions of a continuous 
tangential and a discontinuous normal electric field at the 
border of the cavity r = R: 



<9$ m 
dr 



d<P c 



89 



(A2) 



dr 



By expanding the direct term l/|r — d\ of the internal 
electrostatic field ( & ln at the cavity wall in Legendre poly- 
nomials, solving these boundary conditions is straightfor- 
ward and one finds: 



<J> m = 



<j> out 



Q 

47T£ 

Qk 



1 



\r-d\ 

-^21 + 1 d l 



(1 - 2k) 



l + l d l r 



47re' I + k r 



l+T Pi(cos6) 



R 2l + 



j- Pi (cos ( 



(A3) 



where we introduced k = e'/(e + e'). This result can be 
simplified by rewriting the fraction appearing inside the 
summations: 



<I> m = 



Q_ 

Aire 



1 



d\ 



(l-2/c)(l-«) 
kR 



1 j 



d l r 



1 + kR 21 



Pi (cos t 



(A4) 
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$ out =- 



47re' 



1 - 2k 



Kf 



EK 
TZ 



I + k r l 



Pdcost 



(A5) 



In both cases the first summation represents the expan- 
sion in Legendre polynomials of an inverse distance. In 
the latter this is the distance \r — d\ to the original charge. 
In the former, however, this is the distance to a location 
D == (R/d) 2 d outside the cavity. The second summa- 
tion can be simplified by expanding the Legendre poly- 
nomial in terms of exp(zn6>) (See Ref. 0, Eq.[8. 911.4]). 
This results in: 



1 

47T£ 

0(1 



Q 



\r-d\ 
- 2k 



(1 - 2k 



Q{R/d) 
\r-D\ 



(A6) 



Aire'R 



$ out 



I 2kQ 
Aire' \r-d\ 
Q(l - 2k 



+ 



1 1 



(A7) 



47lV; ^,-,-,1 + ^,2,0 



where we introduced x = r/D and y = d/r, while F\ is a 
hyper-geometric function in two variables (See Ref. [24| . 
Eq. [9. 180.1]). Note that the location D corresponds to 
the position where a single image charge should be placed 
in the case of metallic boundary conditions |l5j . 

Although one can in principle calculate or tabulate 
the hyper-geometric function, one can show that in the 
present case of a confined dipolar fluid this term can 
safely be neglected. In doing so, we approximate the 
induced electrostatic potential field by a single external 
image charge. Note that this approximation is exact in 
the case of a vacuum outside {k = 1/2) and in the case 
of metallic boundary conditions (e' — * oc). 

There are two intuitive arguments for making this ap- 
proximation. Firstly, the field, except near the cavity 



wall, is mainly determined by the local charge distribu- 
tion inside the cavity, rather than the image charge dis- 
tribution arising from the polarisation. The second rea- 
son is that, since we consider a fluid of extended dipolcs 
where the charge separation is roughly a third of the 
particle diameter, the two neglected parts of the induced 
charge distribution of both charges that form a dipole 
will cancel to a large extent. 

In order to illustrate that the error made by the ap- 
proximation is indeed small, we place a unit charge at a 
distance d/R = 0.9 from the origin in a cavity of radius 
R = 2.5 nm and k = 0.8, and measure the reduced force 
F* in terms of the force between two unit charges at a 
distance a. A second charge is put at various distances 

r/R from the origin of the cavity, the result of which is 
o.oi , ,— 



0.008 ■ 



G O r/R = 0.1 

□ □ r/R = 0.3 

O O r/R = 0.5 

A — A r/R = 0.7 
O d r/R = 0.9 



0.006 



< 



0.004 



0.002( f- 




FIG. 17: The absolute error in the reduced force AF* of the 
force of a unit charge at a distance d/R — 0.9 from the origin 
on a unit charge at a position (r, 9) inside a cavity with radius 
R — 2.5 nm and k = 0.8. 



shown in Fig. 1171 The error increases both with decreas- 
ing cavity radius, and on approaching the cavity wall 
with one or both charges. Note, that although the abso- 
lute error increases also when the charges get closer, the 
relative error in that case decreases due to the singular 
behaviour of the direct interaction between the charges. 
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